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For a coronagraph to detect faint exoplanets, it will require focal plane 
wavefront control techniques to continue reaching smaller angular separations 
and higher contrast levels. These correction algorithms are iterative and the 
control methods need an estimate of the electric field at the science camera, 

o , which requires nearly all of the images taken for the correction. The best way 

to make such algorithms the least disruptive to science exposures is to reduce 
the number required to estimate the field. We demonstrate a Kalman filter 
i estimator that uses prior knowledge to create the estimate of the electric 

i field, dramatically reducing the number of exposures required to estimate 

> the image plane electric field while stabilizing the suppression against poor 

CN signal-to- noise (SNR). In addition to a significant reduction in exposures, 

00 

we discuss the relative merit of this algorithm to estimation schemes that 
i~ h do not incorporate prior state estimate history, particularly in regard to 

^ estimate error and covariance. Ultimately the filter will lead to an adap- 

tive algorithm which can estimate physical parameters in the laboratory for 
robustness to variance in the optical train. © 2013 Optical Society of America 

OCIS codes: 110.1080, 350.1260, 100.5070, 110.6770, 350.1270. 



1. Introduction 

The desire to directly image extrasolar terrestrial planets has motivated much research into 
space-based missions. One approach proposed for direct imaging in visible to near-infrared 
light is a coronagraph, which use internal masks and stops to change the point spread function 
of the telescope, creating regions in the image of high contrast where a dim planet can be seen. 
However, coronagraphs possess an extreme sensitivity to wavefront aberrations generated by 
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the errors in the system optics. This necessitates wavefront control algorithms to correct for 
the aberrations and relax manufacturing tolerances and stability requirements within the 
observatory. Of the two, stability is of particular importance as noted by Shaklan [1,2]. In 
this paper we discuss the challenges associated with wavefront estimation and control using 
deformable mirrors (DMs) in a coronagraphic imager, and how we can improve the robustness 
of such a system to observatory instability. Advances in these correction algorithms have 
primarily focused on development of the controller, by choosing some criterion that decides 
how best to suppress aberrations given an estimate of the electric field at that point in time. 
This estimate is found by modulating the system in a predicable fashion so that there is 
adequate phase modulation to solve for the expected value of the electric field at the image 
plane. Whether this modulation is from observations at multiple planes, or by modulating the 
the surface of a DM, the observations are made in the presence of the all prior control signals 
sent to the DMs. Since existing estimation schemes do not predict the effect of a perturbation 
induced by the DM they must abandon all prior knowledge of the electric field estimate, and 
as such require a large number of images to reconstruct the estimate. By utilizing both prior 
estimates and the control history we develop a method that requires fewer images to update 
the estimate, simultaneously improving the efficiency of the correction algorithm and its 
robustness to observatory instability. 

2. Experimental Setup: Princeton High Contrast Imaging Laboratory 



Fig. 1. Optical Layout of the Princeton HCIL. Collimated light is incident 
on two DMs in series, which propagates through a Shaped Pupil, the core of 
the PSF is removed with an image plane mask, and the 90° search areas 3X6 
reimaged on the final camera. 

The High Contrast Imaging Laboratory (HCIL) at Princeton tests coronagraphs and wave- 
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front control algorithms for quasi-static speckle suppression. The collimating optic is a six 
inch off- axis parabola (OAP) followed by two DMs in series and a coronagraph, which is 
imaged with a second six inch OAP (Fig. 1). We use a shaped pupil coronagraph, shown in 
Fig. 2(a), and described in detail in Belikov et al. [3]. This coronagraph produces a discovery 
space with a theoretical contrast of 3.3 x 1CT 10 in two 90° regions as shown in Fig. 2(b). At the 
Princeton HCIL, the aberrations in the system result in an uncorrected average contrast of 
~ 1 x 10~ 4 in the area immediately surrounding the core of the point spread function (PSF), 
which agrees with the simulations shown in Fig. 2(d). Since the coronagraph is a binary 
mask, its contrast performance is fundamentally achromatic, subject only to the physical 
scaling of the PSF with wavelength. 




(mm) 

(a) Shaped Pupil 



(b) Ideal PSF 



\bt;iratt!(l PSF from Ripph;3 




(c) Aberrated Pupil (d) Aberrated PSF 



Fig. 2. (a) A shaped pupil, (b) The ideal PSF from a system using a shaped- 
pupil coronagraph. (c) Shaped Pupil with aberrations generated by Fresnel 
propagating the measured nominal shapes of the DMs to the pupil plane. 
Other sources of aberrations are not included because they have not been 
measured, (d) The PSF of the shaped pupil with the simulated aberrations. 
The figures are in a log scale, and the log of contrast is shown in the colorbars. 
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3. Modeling the HCIL 

The estimation and control algorithms in this paper are dependent on a model of the opti- 
cal system, specifically the transfer function between the DM and image plane. We detail a 
common approach to this mapping, where the field is linearized to describe the problem as 
discrete perturbations beginning at the pupil plane [4-9]. Most commonly, the field is lin- 
earized about the current DM surface shape. In the Princeton HCIL neither DM is conjugate 
to the pupil plane, requiring that we account for this propagation in our model. Beginning 
with arbitrary complex aberrated field, g(u,v), incident on arbitrary aperture, A(u,v), and 
a phase perturbation induced by the DM surface, <fi(u, v), the field at the DM plane is given 
by 

E (u,v) =A(u,v)(l + g(u,v))eWw). (1) 

The DM perturbation commanded by the controller, <fi(u, v), will be added to the prior (also 
referred to as "nominal") DM shape, <f>o(u,v), so that <p = + 0o- Linearizing about the 
nominal DM shape, the field at the pupil plane is given by 

E (u,v) = A(u,v)e i<t,0 (l + g(u,v)+i(l)(u,v)). (2) 

The resultant image plane electric field, E im (x,y), is described by an arbitrary linear 
operator, C{-}, as 

E im (x,y) =C{E (u,v)} (3) 
= C{A(u, v)e^ } + C{A(u, v)e**°g(u, v)} + iC{A(u, v)e uh <j>{u, v)}. (4) 

The intensity distribution of this field is given by 



C{A(u, v)e**°} + C{A(u, v)e^°g(u, v)} + iC{A(u, v)e i<l>0 <f)(u, v)}. 



(5) 



Rewriting Eq. 5 as an inner product and absorbing e l< ^° into C{-}, the intensity distribution 
in the image plane is represented as 

I im =< > +2K{< C{A(1 + g)},iC{A<f>} >} 

+ <C{A(l + g)},C{A(l+g)}>. 

It is important to note that to re-linearize about the current DM shape at any control step, 
k, C{-} must be re-evaluated. Fortunately, this is as simple as convolving C{e l ^ >k - 1 } with each 
component of Eq. 6. 

For the estimator we impose a scalar inner product to each element in the two-dimensional 
plane, (x,y), to compute the intensity distribution across the image plane. For the control 
algorithm we seek an average contrast value inside the dark hole as a metric for the control 
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law. Thus, we impose a matrix inner product to describe the image plane electric field, 
C{Ag}, and the DM control effect, C{A(f)}, our state and control variables respectively, as 
column matrices. These are formed by stacking the columns of the two-dimensional fields 
to create a single column matrix. However, we have yet to parameterize the DM commands 
into a control matrix, u. At the moment we could solve for C{A(f)} but what we seek are 
actuation commands for the DM, not its field at the image plane. To solve for the actuator 
commands, we introduce a physical model of the DM to directly optimize the amplitude 
of each DM actuator. Letting H(x, y) be the height of the DM surface, the resulting phase 
perturbation induced by the DM, <p(u,v), is 

(f>(u,v) = — H(x,y). (7) 

For control, we wish to describe H(x, y) as a combination of the two-dimensional height 
maps imposed by each actuator. Since we are using a DM with a continuous face sheet, the 
contribution of any actuator will be highly localized but still deforms the entire DM surface. 
As a result, we must describe the contribution of the q th actuator as a two-dimensional phase 
map over the entire plane of the DM surface, h q (x,y). The combination of all actuators is 
nonlinear, and very complicated to compute [10], but in the presence of small aberrations 
our deformations will be small. Thus, we assume that superposition of actuators is valid and 
describe H(x,y) as a sum of h q (x,y) over all actuators, N act . The phase contribution of the 
DM is then given by 

2,71 ^ act 

(f)(u,v) = — y2h q (u,v). (8) 

Finally, we wish to make our control matrix, u, a column matrix made up of the control 
signal from each actuator, u q . To do so we describe h q (u,v) as a characteristic shape with 
unitary amplitude, commonly referred to as an influence function, f q (u,v). To find h g (u,v) 
we simply multiply f q (u,v) by the control amplitude, a g . Describing h q (u,v) with influence 
functions, the phase perturbation induced by the DM is 

2tt Nact 

(f)(u, v) = — V a q f q (u, v), (9) 

which sums the q th 2-D phase map, or influence function f q (u,v), for all N act actuators to 
reconstruct (f>(u,v). The strength of each influence function is determined by a q . 

Substituting Eq. 9 into Eq. 6, we write C{A(f)} = C{Af}u, where u is a column matrix 
of actuator strengths, u — [ai . . . ak] T , and C{Af} can be written as a matrix of dimension 
Npi X x N act . To simplify the notation we define this matrix as the control effect matrix, 

G = C{Af} e [N pix x N act ]. (10) 
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This allows us to write 

< C{A<j)}, C{A<f)} >= u T G*Gu, (11) 

where * denotes the conjugate transpose. Applying the matrix form of the control amplitudes 
to Eq. 6, the scalar value for the average intensity in the dark hole, Idh, is given by 

47T 2 m All m 

Idh = —riruMu + — u T %{b] + d. (12) 

A A A 

Where 

M =< C{Af}, C{Af} >= G*G (13) 
b =< C{A(1 + g)}, C{Af} >= G*C{A(1 + g)} (14) 
d =< C{A(1 + g)}, C{A(1 + g)}>= C{A(1 + g)}*C{A(l + g)}. (15) 

Conceptually, d is the column matrix of the intensity contribution from the aberrated field, 
b is a matrix representing the interaction of the DM electric field with the aberrated field, 
and M describes the additive contribution of the DM to the image plane intensity. Having 
represented Idh in a quadratic form we can use Eq. 12 to produce an optimal control strategy, 
but we must first estimate C{A(1+ g)} to compute b at each control step. Doing so optimally 
is the main topic of this paper. 

4. Stroke Minimization in Monochromatic Light 

In general, focal plane wavefront correction methods are broken down into a wavefront esti- 
mation step followed by a control step where the DM surface height is determined. In this 
paper we use the stroke minimizing algorithm described in Pueyo et al. [8] as the controller to 
test the estimation schemes. Stroke minimization corrects the wavefront by minimizing the 
actuator stroke on the DMs subject to a target contrast value [8]. Using a single deformable 
mirror to control the field, these algorithms are capable of correcting both amplitude and 
phase aberrations on a single side of the image plane. Pueyo et al. [8] showed that, to 
first order, two DMs in series are capable of correcting both amplitude and phase aberra- 
tions symmetrically about the optical axis, doubling the discovery space in the image plane. 
Physically, this is due to the phase-to-amplitude mixing from propagating the field between 
non-conjugate planes (the first DM to the second). Expressing the actuator amplitudes of 
both DMs as a vector, u, the optimization problem is written as 



minimize > al = u T u 



N 

ti (16) 



subject to Idh < 10 c ■ 
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Using the quadratic form for I DH in Eq. 12, and applying the constraint to the minimization 
via the Lagrange multiplier, /x, the quadratic cost function is given by 



We find the optimal actuator commands via a line search on \i to minimize the augmented 
cost function, Eq. 17. Since this is a quadratic subprogram of the full nonlinear problem 
we can iterate to reach any target contrast [8]. In addition to regularizing the problem of 
minimizing the contrast in the search area, minimizing the stroke has the added advantage of 
keeping the actuation small and thus within the linear approximation. If the DM model and 
its transformation to the electric field (embedded in the M matrix) were perfectly known, 
the achievable monochromatic contrast would be limited only by estimation error. Note 
that using the derivation shown in §3 we have only explicitly derived the control law for a 
single DM. However, Pueyo et al. [8] have shown that M and b only need to be augmented 
to account for additional DMs at non-conjugate planes. Thus the form of the control law 
remains the same, and only the dimension of the matrices change. 

5. Least-Squares Estimation Using Pairwise Measurements 

To date, almost all high-contrast wavefront control approaches use the DM diversity algo- 
rithm developed by Give'on et al. [6] to estimate the wavefront. Up to this point DM diversity 
had given the best results at the Princeton HCIL in both monochromatic and broadband 
light [8,11], making it the baseline of comparison for the Kalman filter estimation scheme. 
The DM diversity algorithm solves for the electric field by modulating a single DM while 
measuring the variation in intensity at the image plane. We then solve for the field by con- 
structing an observation matrix based on a model of how that DM perturbs the electric field 
at the image plane. 

The linearized interaction of the DM actuation with the aberrated electric field can be 
written in matrix form by taking difference images using j pre-determined shapes with 
amplitudes prescribed by the normalized intensity of the aberrated field [12,13]. The image 
J. + is taken with one deformable mirror shape, 6~, while I~ is the image taken with the 
negative of that shape, —(f>j, applied to the deformable mirror. Applying Eq. 6, the resulting 
measurement is given by 




(17) 



The corresponding optimal command that minimizes Eq. 17 is given by 




(18) 



/+_/-= M{< C{A(1 + g)},iC{A<P} >}. 



(19) 
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Taking j measurements, each with a different conjugate pair of DM shapes, a vector of noisy 
measurements at the k th time step is constructed as 



z k 



L A j,k 



(20) 



We construct the observation matrix, H k , so that it contains the real and imaginary parts of 
the j th DM perturbation, C{A(pj}, in each row. With j DM probe shapes, the observation 
matrix at the k th time step is given by 



^{C{A^ k }} 9f{C{A0 lifc }} 



^{C{A^ k }} Q{C{A<f> jtk }} 



(21) 



Having defined the form of H k , it follows from Eq. 19 that the image plane electric field is 
separated into real and imaginary parts. The state vector describing the image plane electric 
field of a specific pixel at the k th time step is given by 



X{C{Ag k }} 
*{C{Ag k }} 



Assuming our measurement noise, n k , is additive the observation is given by 

z k = H k x k + n k . 



(22) 



(23) 



The true electric field, x k , is the unknown in Eq. 23. Using the observation matrix, H k , we 
seek an estimate of the electric field, x k , based on the noisy intensity measurements defined 
in Eq. 20. The estimator will attempt to minimize the error between the estimated and true 
state. Since the control system does not have access to the true state, the metric defining 
this error is necessarily based off the measurement residual, H k x k — z k . Given a batch of 
measurements at the current time step, k, the cost function used to define the DM diversity 
algorithm is given by 

J = -[H k x k - z k } T R k l [H k x k - z k ], (24) 
where we have weighted the cost function by the measurement covariance, defined as 



R k = E[n k nl\. 



(25) 



The solution that minimizes the weighted quadratic form of the estimated measurement 
residual is given by the left pseudo- inverse of H k , 



Xk 



(H k R k H k ) H k R k z k 



(26) 
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If the measurement noise is assumed to be white, gaussian noise with a standard deviation 
of <7fc, the measurement covariance is given by 



(27) 



and the state estimate provided by calculating the left-pseudo inverse reduces to 



x k = {H k H k ) 1 H k r z k . 



(28) 



With an estimator in hand we wish to quantify its uncertainty by computing its covariance, 
defined as the expected value of the square estimate error, x k — x k . Looking back to the 
weighted form of the estimator in Eq. 26, the covariance of the estimate is given by 



Since this estimator relies on a new batch of measurements at each control step, the covari- 
ance also resets after every control step. Thus the uncertainty in the estimate is tied to the 
SNR and the quality of the observation (the meaning of which is discussed further in §10) 
in that particular set of measurements. 

With this formulation of the estimator, we can estimate the real and imaginary parts of 
the aberrated field at each pixel in the image plane with least-squares minimal error on the 
measurement residual, and by extension minimize the error between the estimated and true 
state of the image plane electric field. A minimum of 2 image pairs must be used to produce 
a unique estimate that minimizes the cost function in Eq. 24. The probe shapes rely on an 
analytical form intended to evenly modulate two rectangular regions of width w x and height 
w y , offset from the optical axis by (a, b) [6,14]. Appealing to the Fourier convolution theorem, 
the DM perturbation that will modulate the image plane in this manner is given by 



where 6 U and 6 V are arbitrary phase shifts that are changed for each probe pair. Practically, 
we find that 4 probe pairs must be used to get a good enough estimate at the Princeton 
HCIL. This is partially a result of the ad-hoc nature with which the phase shift of each 
probe pair may be computed, and each phase shift is not guaranteed to adequately modulate 
the aberrated field. Consequently, 8 images are taken per iteration to estimate the electric 
field when using the DM diversity algorithm. An optimal probing scheme that solves this 
limitation is discussed in §10. 

Experimental results using the DM diversity estimator with the stroke minimization al- 
gorithm are shown in Fig. 3. The laboratory starts at an initial contrast of 1.23 x 10~ 4 , 



P k = E [(x k - x k ){x k - x k ) T ] 



(29) 
(30) 



(p(u,v) 



sinc(w x u) smc(w y v) cos(aw + 6 U ) cos(bv + 9 V ), 



(31) 
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Fig. 3. Experimental results of sequential DM correction using the DM diver- 
sity estimation algorithm. The dark hole is a square opening from 7-10 x -2-2 
X/D on both sides of the image plane, (a) The aberrated image, (b) Contrast 
plot, (c) The corrected image. Image units are log(contrast). 

Fig. 3(a). Using the least-squares estimation technique it is capable of reaching an average 
contrast of 2.3 x 10~ 7 in a (7-10)x(-2-2) X/D region within 30 iterations (Fig. 3(c)) on both 
sides of the image plane, a unique capability that is a result of the two deformable mirrors in 
the system. In 20 iterations of the algorithm, requiring a total of 160 estimation exposures, 
the system reached a contrast level of 3.5 x 10~ 7 . 

6. Constructing the Kalman Filter Estimator 

The DM diversity algorithm [13] is quite effective, but is limited in that it does not close the 
loop on the state estimate, as shown in Fig. 4. Therefore all state estimate information, x, 
acquired about the electric field in the prior control step is lost, requiring a full set of probe 
images to estimate the field at each time step. In addition to being very costly with regard 
to exposures, the measurements will become progressively noisier as higher contrast levels 
are reached. Including feedback of the state estimate at each iteration will add a degree of 
robustness to new, noisy measurements and will reduce the number of exposures required 
to update the electric field estimate at each control step. Thus, incorporating prior history 
means the estimator will reach its final contrast target with a minimal set of exposures and 
will be less sensitive to poor SNR measurements as the controller suppresses the field. To 
include state feedback we must extrapolate the controller's perturbation to field from the 
last control step using the model in §3. Since the estimate will now be formed by combining 
an extrapolation of the prior estimate with a smaller set of measurement updates, we must 
consider the relative effect of process and detector noise to optimally combine the two. This 
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Fig. 4. Block diagram of a standard FPWC control loop. At time step k, only 
the intensity measurements, Zk, provide any feedback to estimate the current 
state, Xk- The dashed lines show additional feedback included in this paper 
from the prior electric field (or state) estimate, Xk, and the control signal, Uk- 



is exactly the problem a discrete time Kalman filter solves. The noisy measurement updates, 
Zk — Uk + n k , will still be difference images of probe pairs with a covariance defined by Eq. 25. 
The conjugate pairs allow us to construct a linear observation matrix, H k . If we were not in 
a low aberration regime our observer would have to be nonlinear. This is not impossible for 
a Kalman filter, but can make it highly biased [15] and computationally expensive. 

Like the DM diversity estimator [6,13], the Kalman filter produces an estimate with least- 
squares minimal error. Since the Kalman filter operates on the estimate in closed loop, the 
weighted cost function given in Eq. 24 will not adequately represent the error contributions 
in the system. To account for the noise included in forward propagating our prior state 
estimate history in time, we must also include an estimate of the state covariance, defined in 
Eq. 29. Since the covariance as yet has not been updated with a measurement, we denote it 
as Pfc(— ), the minus sign indicating that it is an extrapolation from our model. Defining the 
error as both the difference between the noisy observation and the estimated observation, 
HkXk{+)—z, and the difference between our current estimate and the estimate extrapolation, 
%k ~ Xk{—), the new quadratic cost function is given by 



J = [Xk ~ X k {-)] T Pk 



[Xk - x k { 



[H k x k - z k ] T R k 1 [H k x k 



Zk 



(32) 



We formulate the cost in matrix form as 



J 



1 


x k - x k {-) 


T 


~P k {-) o " 


-l 


x k - x k {-) 


2 


H k x k — z k 




R k _ 




HkX k — z k 



HkX k — z k ) T R k 1 (H k x k — Zk), 



(33) 
(34) 
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where we have now defined a new set of augmented matrices as 



Hh 



Rk 



x 

x k (-) 
z k 

'Pk(-) 
R k 



(35) 
(36) 
(37) 



Evaluating the partial derivative of Eq. 34 the state estimate update is given by 

x k {+) = x k (-) + P k {-)Hl [H k P k {-)Hl + R k ]~ l [z k - H k x k (-)} . (38) 

From Eq. 38, we define the optimal gain to be 

K k = P k {-)H T k [H k P k {-)Hl + R k }~\ (39) 

Finally, we update the state covariance estimate, P k (+), by applying Eq. 38 to the expected 
value of the covariance, 



P k (+) = E[(x k (+) - x k )(x k (+) - x k ) 



(40) 



which gives 

Pk(+) = [Pk(-)- 1 + HlR- k 'H k \~\ (41) 

To complete the filter we need to propagate the prior estimate, x k -i(+), to the current 
time step. The filter extrapolates to the current state estimate, x k (—), by applying a time 
update to the prior state estimate via the state transition matrix, <& k -i, and numerically 
propagating the control output from stroke minimization at the prior iteration, u k -i, via a 
linear transformation described by r fe _i. We also have a disturbance from the process noise, 
u>jfc_i, which is propagated to the current state of the electric field via the linear transfor- 
mation, A fe _!. Assuming these components are additive, the state estimate extrapolation is 
given by 

x k (-) = $jfc-ix fe _i(+) + r fc _iw fc _i + A fc _iiUfc_i. (42) 

To determine these matrices, we apply the linearized optical model developed in §3. Using 
a linearized model avoids generating arbitrary bias in the estimate at each pixel, a common 
problem with a nonlinear filter [16]. The first step in propagating the state forward in time 
is to update any dynamic variation between the discrete time steps with the state transition 
matrix, In this system, captures any variation of the field due to temperature 

fluctuations, vibration, or air turbulence that perturb the optical system. To simplify the 
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model, we recognize that there is no reliable way to measure or approximate small changes 
in the optical system over time with alternate sensors; we assume that the state remains 
constant between control steps, making the state transition matrix the identity, $ k -i = 
$ = X. However, it is very important to note that $fc_i provides a simple mechanism to 
incorporate time updates from alternative high speed sensors, such as residuals measured by 
the adaptive optics wavefront sensor, between image plane measurements. Each submatrix for 
T, shown in Table 1, is of dimension 2 x 2Nom- Every pair of rows is found by separating the 
real and imaginary parts of Eq. 10 into individual rows. Making the standard assumption 
that the process noise is gaussian white noise, the expected value of the state when we 
extrapolate is 

x fc (-) = $ fc _ix fc _i(+) + r fc _iit fc _i. (43) 
It's associated covariance extrapolation is then given by 

p fc (-) = Sfc-iiW+^Li + Qk-i- (44) 

Combining Eq. 38, Eq. 39, and, Eq. 41 with the extrapolation equations, this form of the 
filter consists of five equations that describe the state estimate extrapolation, covariance esti- 
mate extrapolation, filter gain computation, state estimate update, and covariance estimate 
update at the k th iteration [15]: 



x fc (-) = $jfc_ix fe _i(+) + r fc _itt fc _i. (45) 

PkH = + Qk-i (46) 

K k = P k {-)Hl [H k P k {-)Hl + R k ] 1 (47) 

x k {+) = x k {-) + K k [z k - H k x k {-)] (48) 

P k (+) = [Pki-r 1 + H T k R k l H k y l (49) 



The focal plane measurements z k are identical to that of §3, and are constructed into a 
column vector of j difference images for n pixels. Likewise H k takes on a similar form, and 
is a matrix constructed from the effect of a specific deformable mirror shape <f)j tk on the 
real and imaginary parts of the electric field in the image plane. Finally, we compute the 
covariance update, P k (+), based on the added noise from the new measurements. As with 
the DM diversity estimator, the estimated state is a column vector of the real and imaginary 
parts of the electric field at each pixel of the dark hole. The control signal u remains a 
column vector of the DM actuator strengths as defined in §3, with DM1 commands being 
stacked on top of the DM2 commands. The dimension and form of the rest of the filter 
equations follows from these. They are all listed in detail in Table 1 and Table 2. Since we 
are only considering process noise at the DMs, the process disturbance w is a vertical stack 
of the variance expected from each actuator. The only remaining step critical to the filter's 
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performance is initializing the covariance, P . In our system this cannot be measured, so we 
must initialize with a reasonable guess. 

The Kalman filter gain, Eq. 47, optimally combines the prior estimate history with 
measurement updates to minimize the total error contributions based on the expected state 
and measurement covariance. A fundamental property of the Kalman filter is that the opti- 
mal gain, Eq. 47, is not based on measurements, but rather estimates of the state covariance, 
Pk(-), process noise from the actuation Qk-i, and sensor noise Rk- This means that the op- 
timally of the estimate is closely related to the accuracy and form of these matrices; this 
will be discussed in §8. Much like the batch process method the purpose of the gain in 
the Kalman filter is to minimize a quadratic cost function, Eq. 32, but it is also subject to 
the constraining dynamic equations defining Xk(-) and Pk{— )■ However, looking at Eq. 35 
there is a major advantage of the Kalman filter in its minimization of the cost function. For 
the batch process method, the observation matrix will be underdetermined if less than 2 
probe pairs are applied. Thus the solution to cost function in Eq. 24 is non-unique. However 
the solution for the Kalman filter is guaranteed to be critically determined even with no 
measurement update, and Hk only requires a single measurement to provide an overdeter- 
mined solution to the cost function in Eq. 34. Thus, the Kalman filter is formulated in such a 
way that it solves a least squares error, left pseudo-inverse problem, regardless of the number 
of measurements taken. This gives us the freedom to minimize the total number of exposures 
required to estimate the field with enough accuracy to reach the target contrast level. 

Referring back to Eq. 30, we see that in both cases the covariance depends on the current 
observation, but the Kalman filter is also a function of the prior state covariance. Looking 
at Eq. 49 the final term, H]:R^ l H k) is guaranteed to be positive definite, meaning that 
additional measurements can do nothing but reduce the magnitude of the covariance, Pk(+)- 
In the event of a measurement with a poor SNR the covariance may not get better, but it 
is guaranteed not to get worse because H^R^Hk is positive definite. This is indicative of 
the estimator's robustness to bad measurements at later time steps. In contrast, Eq. 30 
demonstrates that if the left-pseudo inverse solution is given measurements with a poor SNR 
the new estimate will have a large covariance at this time step. In this case the control 
will not be effective, which is why we often see jumps in contrast when the DM Diversity 
estimator reaches low contrast levels, as seen in Fig. 3. In the case of the Kalman filter, the 
high covariance contributed by a poor measurement update is dampened by the contribution 
of prior covariance estimates via Pk(-), stabilizing the state estimate and its covariance in 
the event of a bad measurement. Since we cannot guarantee that a probe will provide a 
high enough SNR, particularly at high contrast levels where we have reduced the number of 
photons per speckle, this is an extremely attractive component of the Kalman filter estimator. 
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Matrix 



Dimension 



$ =X 



(2 • N p i xe i s ) x (2 • Npi xe i s ) 



r = 



^■{Gdmi} ^{GdM2} 
^{Gdmi} ^{GdM2} 

^■{Gdmi} ^{GdM2} 



/ 



Hk = diag 



^{G DM2(f>l,k} ^{G DM2<f>l,k} 



1 \ 



\^ lU{G DM2 (f)j,k} ^{G D M24>j,k}\ J 

K k is computed 



u = 



DM1 
DM2 



z = 



it - ir 



'it - 1 1 



(2 • N pixels ) x (2 • N DM ) 



{■^pixels) ' {-^pairs) X (2 • -^Vpixe/s) 



(2 ' Npixeis) X (.-^pixels) ' (,-Npairs) 



(2 • Ndm) x 1 



\Npixels) ' \Npairs) X 1 



X = 



(2 • iVp i;Eeis ) x 1 



Table 1. Definition of all propagation Matrices. Ndm is the number of actuators 
on a single DM, N pixe i s is the number of pixels in the area targeted for dark 
hole generation, and N pa i rs is the number of image pairs taken while applying 
positive and negative shapes to the deformable mirror, j indexes the probe 
shape, k the discrete time step, and n the pixel the dark hole. 
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Matrix 



Dimension 



P = E[(x Q - x )(x - x ) T ] (2 • N pixels ) x (2 • Npi xels ) 

Q k = A£K<]A T (2 • N pixels ) x (2 • N^,) 

A = r (2 • iV pi;Eeis ) x (2 • 7V DM ) 



&DM1 
&DM2 



(2-N, 



DM, 



X 1 



Table 2. Definition of all noise Matrices. N DM is the number of actuators on 
a single DM, N pi xeis is the number of pixels in the area targeted for dark 
hole generation, and N pairs is the number of image pairs taken while applying 
positive and negative shapes to the deformable mirror, j indexes the probe 
shape, k the discrete time step, and n the pixel the dark hole. 



7. Iterative Kalman Filter 

Since we had to linearize the field to obtain Tk~i and Hk, we iterate the filter on itself 
to account for some of the nonlinearity not captured in the initial linearization. We do so 
by feeding the newly computed state and covariance update Pk(+) back into the 

filter again, setting u^-\ to zero. Applying feedback to the filter in this manner is purely 
computational, and comes at no cost with regard to exposures. For sufficiently small control 
this additional computation will account for nonlinearity in the actuation and better filter 
noise in the system, limited only by the accuracy of the observation matrix, H k . Improving 
the accuracy at each time step means the correction algorithm will require fewer iterations, 
and hence fewer exposures, to reach its final contrast target. For a severely photon-limited 
exoplanet imager, such a reduction minimizes the time required to suppress the speckle field. 
Following a notation similar to Gelb [16], the p th iteration (for p > 2) of numerical feedback 
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into the iterative Kalman filter at the k th control step is 



%p,k{ — ) — &p-l,k(+) 

Pp,k{~) = ®p-l,kPp-l,k(+)%-l,k + Qp-l,k 



(50) 
(51) 
(52) 
(53) 
(54) 




-i 



The power of iterating the filter lies in what we are fundamentally trying to achieve. For a 
successful control signal, we will have suppressed the field. This means that the magnitude 
of the probe signal will be lower than the control perturbation. This guarantees that Hk 
will better satisfy the linearity condition than Tu. As a result, if we iterate the filter on 
itself during a given control step we can use the discrepancy between the image predicted by 
Hk%k{ J r) and the measurements, Zk, in Eq. 53 to filter out any error due to nonlinear terms 
not accounted for in T. In this way, we can accommodate a small amount of nonlinearity in 
our extrapolation of the state without having to resort to a nonlinear, or extended, Kalman 
filter. This means that in the work presented here we did not have to re-linearize within a 
single control step, k, as would be the case for an iterative extended Kalman filter (IEKF). 
It also avoided having to concern ourselves with any bias introduced into the estimate by 
a nonlinear filter. While we have chosen not to in this case, we can move to an IEKF in 
the future by simply re-linearizing about the current DM shape and iterating on p until 
the estimate converges. Something very important to point out is that unlike many IEKF 
problems, this has a unique flexibility with regard to the linearization. A conventional IEKF 
has been linearized about the state, Xk, meaning that the new matrices in the IEKF must 
remain within an allowable radius of convergence, oftentimes requiring that the estimate be 
re-linearized at every iteration on p [16]. Looking back to Eq. 2 however, we recall that we 
have actually linearized about the control signal. Since we did not linearize about the state 
and the problem assumes no dynamics at this point, we need only re-linearize IV- 1 and Hk- 
This re-linearization need not wait for the estimate to begin. Once the control at (k — 1) 
is applied we can re-linearize T k _i and H k in parallel with taking new measurements for 
the current time step, z^.Thus for an adequately long exposure, re-linearizing these matrices 
can be scheduled in such a way that they have negligible impact on the time required to 
estimate the field. Additionally, Eq. 46 and Eq. 47 are measurement independent. Thus, their 
computation can be scheduled in parallel to taking the new measurements, further improving 
the efficiency of the estimator with regard to time. 
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8. Sensor and Process Noise 

Two important design parameters for the performance of the filter are the process noise, Qk-i, 
and the sensor noise, Rk. For the laboratory experiments we make reasonable assumptions 
by appealing to physical scaling of the two largest known sources of error in the system. 
The sensor noise will be determined by the dark current and read noise inherent to our 
detector. In a space observatory, and possibly even on the ground, the sensor noise will 
also be determined by photon noise. We assume that the process noise is dominated by 
errors in the commanded DM shapes. Thus, the process noise is isolated to poor knowledge 
of the DM surface, which comes from the inherent nonlinearity in the voltage-to-actuation 
gain as a function of voltage, the variance in this gain from actuator to actuator, and the 
accuracy of the superposition model used to construct the mirror surface that covers the 
32x32 actuator array of the Boston Micromachines kilo-DM. Treating actuation errors as 
additive process noise, the Kalman filter can account for them in a statistical fashion, rather 
than deterministically in a physical model. Since there is no physical reasoning to justify 
varying Q k at each iteration, it will be kept constant throughout the entire control history 
(Qk-i = Q = constant). Two versions of process noise can be considered in this case. The 
first is where there is no correlation between actuators, giving a purely diagonal matrix 
with a magnitude corresponding to the square of the actuation variance, a u . The second 
version of Q has symmetric off-diagonal elements, which treats uncertainty due to inter- 
actuator coupling and errors in the superposition model statistically. As a first step, we will 
not consider inter-actuator coupling to help avoid a poorly conditioned matrix. This helps 
guarantee that the Kalman filter itself will be well behaved. Thus the process noise for the 
filter will be 



Q = alTXT T . (55) 

Following Howell [17] the noise from both the incident light and dark noise in a CCD 
detector follows a Poisson distribution. Since each measurement is a difference of pairwise 
images the noise is zero-mean and will become more Gaussian as the exposure time increases. 
We simplify the noise statistics by assuming it is uncorrelated and constant from pixel to 
pixel. The measurement error covariance, Rk, is thus a diagonal matrix of the mean pixel 
covariance, oqcd, given by 

2 

R = a CCD j ( 5 gN 

/* Npairs^- N pairs ' \ / 

oo 

where I 00 is the peak count rate of the PSF's core (allowing us to describe R in units of con- 
trast). Having appealed to physical scaling in the HCIL, we now have close approximations 
of the true process and sensor noise exhibited in the experiment. 
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9. Experimental Results for the Kalman Filter Estimator in Monochromatic 
Light 

We corrected the field using the Kalman filter estimator using four, three, two, and one pair 
of images as a measurement update to assess the degradation in performance as information 
is lost. We began with four measurements (four image pairs), to compare its performance 
using the same number of measurements as were used in the DM diversity estimator. Using 4 
pairs, the filter achieved a contrast of 4.0 x 10 in (7-10)x(-2-2) X/D symmetric dark holes 
within 20 iterations of the controller, shown in Fig. 5. Note that this used a total of 160 
estimation images, which is the same amount of information available to the DM diversity 
estimator when it achieved a contrast of 3.5 x 10 -7 in 20 iterations. 



Initial Image Exposure Number Final Image: 4 Pairs 




-10 -5 5 10 o 10 20 -10 -5 5 10 

X /D Iteration A /D 



(a) (b) (c) 

Fig. 5. Experimental results of sequential DM correction using the discrete 
time extended Kalman filter with 4 image pairs to build the image plane 
measurement, Zk- The dark hole is a square opening from 7-10 x -2-2 X/D 
on both sides of the image plane, (a) The aberrated image, (b) Contrast plot, 
(c) The corrected image. Image units are log(contrast). 

Reducing the number of image pairs to three, the correction algorithm reached a contrast 
level of 5.0 x 10~ 7 using only 120 estimation images, as shown in Fig.6. After fine tuning 
the covariance and noise matrices, the correction algorithm achieves a contrast of 2.3 x 10~ 7 
in 30 iterations using only two image pairs per iteration, shown in Fig. 7. The results are 
better than the three and four pair cases because we improved the covariance initialization 
and increased the number of numerical iterations (p from §7) of the Kalman filter for a 
single control step. As discussed in §7, numerically iterating the filter is critical to its per- 
formance, despite the fact that no additional measurements are taken, because it accounts 
for nonlinearity not captured when numerically propagating the control signal via IV_i. 
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Fig. 6. Experimental results of sequential DM correction using the discrete 
time extended Kalman filter with 3 image pairs to build the image plane 
measurement, Zk- The dark hole is a square opening from 7-10 x -2-2 X/D 
on both sides of the image plane, (a) The aberrated image, (b) Contrast plot, 
(c) The corrected image. Image units are log(contrast). 




Fig. 7. Experimental results of sequential DM correction using the discrete 
time extended Kalman filter with 2 image pairs to build the image plane 
measurement, Zk- The dark hole is a square opening from 7-10 x -2-2 X/D 
on both sides of the image plane, (a) The aberrated image, (b) Contrast plot, 
(c) The corrected image. Image units are log(contrast). 



Reducing the number of measurements to a single pair we find a very interesting result. The 
quality of the measurement at any particular time step of the algorithm is now dependent on 
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the quality of that particular probe shape. If a probe does not happen to modulate the field 
well Hk drops rank, adding noise to the gain computation and estimate update of Eq. 47 
and Eq. 48. As a result, the controller can actually degrade the contrast in the dark hole at 
that time step. Additionally, a single probe may not modulate a specific location of the field 
well, so we must choose a different probe shape at each iteration to guarantee good coverage 
of the entire dark hole in closed loop. Starting from an aberrated field with an average 
contrast of 9.418 x 10™ 5 , Fig. 8(a), we achieved a contrast of 3.1 x 1CT 7 in 30 iterations and 
2.5 x 10~ 7 in 43 iterations of control, Fig. 8(c). Looking at the contrast plot in Fig. 8(b), 
the sensitivity of a single measurement update to the quality of the probe is very clear. 
What is interesting, however, is that the modulation in contrast damps out over the control 
history. While we do not suppress as quickly with regard to iteration we achieve our ultimate 
contrast levels with fewer total measurements, which is the ultimate performance metric for 
efficiency of the estimation and control algorithm. Thus, even with one measurement update 
per iteration the prior history stabilizes the estimate in the presence of the measurement 
update's poor signal-to-noise at high contrast levels. More encouraging is that these results 
use the analytical probe shapes described by Eq. 31. The choice in phase shifts, 9 U and 9 V , 
is somewhat arbitrary, with no criteria to evaluate how effectively each modulates the field. 
In fact, we had to carefully choose our phase shifts to improve the convergence rate of the 
single probe pair results. If we could choose probe shapes that are guaranteed to modulate 
the dark hole well, we would see a dramatic improvement in the rate of convergence for a 
single measurement update. 

A very promising aspect of this estimation scheme is that its performance did not degrade 
as the amount of measurement data was reduced. With only 86 estimation images it was 
capable of reaching the same final contrast achieved by the DM diversity algorithm in §3, 
which required 240 images to maintain an estimate of the entire control history. Thus by 
making the estimation method more dependent on a model we were able to reduce our need 
to measure deterministic perturbations in the image plane electric field. 

10. Optimal Probes: Using the Control Signal 

The results of §9 use the analytical probe shapes described by Eq. 31. As discussed in §9, 
the choice of the probe shape used for estimation is critical to create a well-posed problem, 
and has been somewhat of an art in the high contrast imaging community. These shapes 
are intended to modulate the field as uniformly as possible so that Hk is full rank. However, 
no formalism has been proposed to determine the "best" shapes to probe the dark hole. In 
any dark hole there are discrete aberrations that are much brighter than others, requiring 
that we apply more amplitude to those spatial frequencies. Conversely the bright speckles 
raise the amplitude of the probe shape, making it too bright to measure the dim speckles 
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Fig. 8. Experimental results of sequential DM correction using the discrete 
time extended Kalman filter with one image pair to build the image plane 
measurement, 2%. The dark hole is a square opening from 7-10 x -2-2 X/D 
on both sides of the image plane, (a) The aberrated image, (b) Contrast plot, 
(c) The corrected image. Image units are log(contrast). 

with an adequate SNR. Excluding this issue, we also cannot truly generate the analytical 
functions, as the actual surface is a sum of individual influence functions. Even the DM 
with the highest actuator density available, the Boston Micromachines 4K-DM©, can only 
approximate each function with 64 actuators. We account for the true shape in our model 
but this shape does not truly probe each pixel in the dark hole with equal weight, which was 
the primary advantage of using the analytical function for a probe shape in the first place. 
Fortunately, we can once again appeal to our mathematical model for estimation and control 
to help determine an adequate probe shape. In closed loop, the control law incorporates the 
influence function model of the surface to determine a shape that necessarily modulates the 
aberrated field in the dark hole. In principle, if we apply the conjugate of the control shape 
we will increase the energy of the aberrated field. Instead of applying probes in addition to 
the control shape, we use the control shape itself (and its negative) to probe the elecric field 
in the dark hole. Thus, we can rely on the controller to compute optimal probe shapes that 
will inherently modulate brighter speckles more strongly than dimmer speckles and we have 
confidence that the shape will adequately perturb the dark hole field. 

Fig. 9 shows preliminary results using the control and its conjugate shape as the probe 
pair for estimating the field with the Kalman filter. The current contrast level is limited to 
2.30 x 10~ 6 on both sides of the image plane. The asymmetry of the dark holes actually hints 
at the reason for not achieving a higher contrast level. In a 2-DM system, the idea that the 
optimal control signal can be used as the probe shape requires that we use both mirrors to 
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Fig. 9. Experimental results of sequential DM correction using the discrete 
time extended Kalman filter with the control effect and it's conjugate shape 
used as the only probe pair for the measurement update, Zk- The dark hole is 
a square opening from 7-10 x -2-2 X/D on both sides of the image plane, (a) 
The aberrated image, (b) Contrast plot, (c) The corrected image. Image units 
are log(contrast). 



probe the field. Since the estimator is currently written assuming one mirror is probing the 
field, we were forced to collapse both control shapes onto the same DM. This means that 
we are not fully perturbing the field as we would have expected, leaving a particular set 
of aberrations unprobed. This is not a problem for a single sided dark hole using one DM, 
where the control and probe surfaces are one and the same. For a two DM system, we simply 
need to reformulate the estimator as a function of both DMs for this concept to be effective. 

11. Conclusions 

We have demonstrated a discrete time iterative Kalman filter to estimate the image plane 
electric field in closed loop. The experimental results only required a single linearization for 
the estimator, and we have mathematically demonstrated a suitable extension to an IEKF 
in the event that we must re-linearize about a new DM shape. This type of progress is 
critical for improving the efficiency and robustness of future coronagraphic missions, thereby 
maximizing the likelihood of planetary detections. We demonstrate the fastest suppression 
to date in the Princeton HCIL, only requiring a single measurement at each iteration. This 
represents a ~70% reduction in the number of images required to estimate the image plane 
electric field. Having increased the speed of the estimator, this makes it a feasible focal plane 
estimation technique for ground-based coronagraphic instruments. The closed loop nature of 
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the estimator is attractive because it also stabilizes the estimate to measurement noise. So 
long as the probe adequately modulates the intended dark hole area, a measurement update 
with poor signal-to-noise does not adversely affect the covariance of the state estimate. 
Having controlled the field with only a single measurement update, we have demonstrated 
that not all probe shapes are best for estimation, motivating us to attempt using the control 
shape as the probing function. Preliminary results for using the control shape as a probe 
signal are promising, and may further reduce the number of required exposures to one image 
per iteration. The Kalman filter also opens up the possibility of adaptive control techniques 
to learn laboratory physical parameters, sensor fusion concepts for integration with extreme 
adaptive optics systems, and bias estimation to gain certainty in planetary detection using 
only the control history. 
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